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in the early 1960's, much effort has been devoted to the 
development of the method itself, while only recently has 
there been any research directed at minimizing the discreti- 
zation error by a proper selection of the element grid. This 
paper examines some recently proposed grid optimization 
techniques and applies them to some one- and two-dimensional 
linear self-adjoint boundary value problems. Guidelines 
reguiring minimal scftwar2 modification are recommended to 
asSist the analyst in obtaining improved finite element 
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I. INTRODUCTION 
The critical concern in any finite element analysis is 
the adequacy of the selected mesh to provide reliabie solu- 
tion results within seme reasonabl2 cost range. The goal of 
finite element grid optimization shen becomes one of 
obtaining maximum solution accuracy for a prescribed anal- 


Mees cost. While this objective is generally not realized in 
today's widesoread use of finite element analysis, the effri- 
Grent computation of solutions with optimal 

become paramount to the engineer as finite els 
mecoeapplz.ed tO increasingly difficult dynamic, nonlinear, 


and evolution problems. 


A. HISTORICAL BACKGROUND 


In the early 1960's, with the help of the high speed 


jer, 


digital computer, finite element methods revolutionize 
problem solving in engineering. Since that time the major 
TeseCarch efforts have concentrated on axpanding the theoret- 


ical basis of the method and extending its application in a 


variety cf fields. Only recently has there been significant 
MeeentlOon directed at minimizing finite element sclution 
a Dts 


smeecs Dy properly defining the e 
@eee a2stributing the nodes a 2 
ensure scme degree of confidence in the solution accuracy 
wer2 poredominantly dependent upon the an 

Judgement and experience, since thare w 


procedures for accomplishing this 
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Mee=emots tcwartrdsS grid optimization h 
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ed with the advent of automatic mesh generators, whic! 
Ss 
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it } 
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v drastically reduced preprocessing SCice 





aeconpmiching tittle ian improving solution accuracy. MTne 
programs automatically construct the element grid, usual 
in auniform manner, after the user merely defines the 
problem and specifies the number of elements to be used. To 
establish convergence and achieve the desired solution accu- 
racy, the user simply repeats th2 analysis using a finer 
mesh of uniformly distributed elements while monitcring such 
convergence indicators as successive solution values at 
common nodal points or the asymptotic increas? in the energy 
content of the mesh. The often disastrous flaw in such a 
practice is that for nearly degenerate problems which 
exhibit very large gradients, the asymptotic range is only 
entered fer an extremely large number of degrees-of-frsedon, 
often beyond computer limitations [Ref. 1]. Pee nase Cace, 
uniform mesh refinement may produce apparent convergence, 
when in fact the solution accuracy is poor. Therefors, the 
need for a finite element grid sptimization precedure is 
clearly evident. 

Phemtirst £ormal attempts at Linite element grid optini- 
meson did not begin unti Snemoaiiy 1970 Vs: These early 
approaches involved the inclusion of the nodal coordina<es 
as dependent variables in the mininization of the potential 
Siecagy £unctional [Ref. 2}. Unfortunately, the resulting 
System of equations is highly nonlinear andthe computa- 
meonal effort involved in its solution is so great that 
Eilat accuracy can be obtained at 4 fraction of the 
expense, Simply by employing a very fine mesh. Clearly, this 
method does not appreach the finite alement grid optimiza- 
meet goal OF achieving a specified solution accuracy for a 
Minimum analysis cost. For this reason, virtually all of the 
gtid optimization technigues since developed are based ona 
"mear-optinum" strategy whereby nearly-optinal solution 
Besults are obtained without the computational inefficiercy 


cE a numerical Optimization analysis. The growing emphasis 
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noSeeenmennadapelvViey, @aprocedure for efficient ceonstruc- 
tion of near->ptimum grids by the iterative application of 
some criterion, based on data already computed from tne 
soluction for a previous grid. This procedure is far more 


efficient than the conventional approach of repeating the 


analysis using successively finer uniform GrLe si. 
Experimental self-adaptive finite element ccdes have 
recently been developed. tarting from the user's initial 


idealization, these programs automatically generate a near- 


optimum grid and solve the resulting equations. 


Be INVESTIGATIVE APPROACH 


In undertaking any numerical optimization task, the 
analyst must first define the objective along with any 
constraints t>9 be imposed upon the objective variables; 
finally a numerical algorithm must be prescribed to perio 
mars optimization, preferably one which wili do so effi- 
ciently for the particular problem. Since the tarm "“optimun" 
most often refers to a solution obtained by nathematical 
Becgbamming, Which is very inerfislent for grid optimiza- 
een, a neat-optimum grid obtained by application of an 
adaptive procedure, henceforth will be termed an optinun 
grid. However, before such a grid can be determined to 
Satisfy the stated cbhjective of obtaining maximum solution 
Meemuracy =2Or a prescribed ccest, che, temms "Accuracy and 
Meest” must be defined; but, more importantly, the eptimiza- 
*ion goals must be specified. This is critical because aqrid 
optimization can be implemeneed in various forms depending 
upon the optimization goals, which wiil, in general, bs 
determined by the original purpose for performing the finite 
element analysis {Ref. 3}. Hon ckanoine, = Lf the purpess. of 
the analysis is to evaluate a local quantity, Such as the 


Maximum value of the denvendent variable or one of it 


WW 





derivatives, then the nodal distribution should be densest 
in the region where that maximum is determined. Eee On toe 
other hand, the interest is on an integral guantizy of the 
dependent variable cver a region of the domain, then the 
nodes should be assigned to achieve a uniform distribution 
of the error over that region. Fon cne welaeose Of ~chis 
investigation, attention wili be focused on *he three finite 
element resultants with the nost Significance in eubid 


soli 
mechanics and other fields in which finite elament analysi 
is employed: the maximum value of the dependent variable, or 
Semsrta2cn: che maxamum Value of the gradient of the soluticn 
and the integral of the solution ovér the domain. 

Mmimords: | tO Tdefinerthe solution accuracy, it wiil be 
necessary tc compare the error in the solution cresultant 
Obtained using an Optimal grid to tha error obtained using 
some baseline grid with the same number cf degrees-cf- 
freedcm. For convenience, the referance grid chosen will be 
a uniform grid, or one with all elements of the Same crder 
and approximately the same size, wita the understanding chat 
no knowledgeable analyst would attenpt to use such a grid in 
the solution of a prcblem with large gradients. 

The definition of anaiysis cost will be greatly simpli- 
fied by making the assumption that it is directly proper- 
e2onal < the number of deqrees-of-fireedom used ‘tc obtain 
mime Solution. Pieneeleaayeene COs. CependsS on many factors, 
sone of man dae are very Cigstac uit =0 quantify. 
Understandabiy, the number of degraes-of-rrseedonm is not the 
sole measure st computational costs, but it appears to be an 
adequate measur of preprocessing and postprocessing costs 
Mermen Cl Teh account for the Major portion of the total 


analysis. 


a 


fem rest idarticn vill concent: 
element grid optimizgation methods fer scl 


v 
Seeucc“rtal mechanics. @hiie this area has doni 





literature on the subject, 
extend equally as well to 


orinciples apply. 


There are two key questions 


~o the 
epeimization: 


adaptive 


application of 


the techniques presented herein 


any E2eld f£6L M@ehich vatriaticnal 


which must be answered prior 


Finite element grid 


(1) What criterion, based on the resuits of an initial 
finite element analysis, should be used to indicate 
re _ s where the original grid is inadequate ? 


(2) What method of grid refinement should be empleyed ? 


Considerable atten 
PeHe next two chapters. 


C. OBJECTIVES 


The objectives of this 


(1) To examine some recently deveioped gri 


techniques which have reached 


appl veat LO. 


(2) To implament some of 


ion of some cne- 


tion will be devoted to 


and two-dinensional 


these dquestions in 


investigation are: 


"N= Seate Of Diace3 eal 


the solu- 
self- 


these techniques in 


linear 


adjoint boundary-value problems. 


(3) To draw a comparison 


ain serms of solution 


among th 


2sé applied techniques 


Sccugacy, analyszs cost, and 


2ase of implementation. 


(4) To reccnmend some 


in obtaining oD 


SELOvIng currently 


so fi tWars. 


“he analyst 
eélemenz solutions 


or easily amendable 
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II. CRITERIA FOR GRID REFINEMENT 

PHeweraMarvy  thecseticai Comecerh in the application of 
adaptive grid optimization is the selection of the refine- 
ment criterion. In cther words, on2 must decide upon which 
Pommcion patraheter, obtained from an initial idealization, 
May most appropriately be used to give some indication as to 
where the initial grid is inadequate and thus needs refine- 
ment. There are several competing proposals concerning the 
most appropriate choice of a refinement criterion. i 
reality, the decision must be based upon such factors as the 
type of problem being solved, which criterion is most prac- 
tically implemented, and whether accuracy is ds¢ssired 
locally, globally, pcintwise, or with respect to an integral 
quantity. The following are some of the more practical 


refinement criteria used in grid optimization at present. 


Ae SOLUTION PARAMETER VARI ATION 


The most Beect, computationally inexpensive, and 
earliest proposed indication of wnere an element grid 
requires refinement is a measure of the variation of some 
solution parameter over the domain. It is only logical chat 
a piecewise polynomial approximation would experience the 
most difficulty in modeling the desired response in «hose 
tegions where the solution or its resultants were either not 
smooth or were characterized by large gradients. Therefore, 
me Dasis cf this criterion is to refine the mesh in those 
areas where a solution parameter varies rapidly, with the 
implication that the optimum mesh is one for which the solu- 
tion parameter variation over e2ach element is uniforn 


throughout the domain. Wics@-@steecOnsedesation ih ‘the 
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pppeeeeesiOneer such @ criterion 15 to tind a scheme for 
Sectcibuting the nodes ¢*5 achieve such a condition. For 
one-dimensional probiems the task is trivial, but one way to 
ensure equal variation over each element in higher dimen- 
Sions is to position the nodes along equidistant contours of 
the chosen parameter. This is precisely the procedures 
prescribed for a practical optimization technigue known as 
contouring. The other consideration is the determination of 
which solution parameter is to be used. MAmtact., Several 
solution parameters have been found +0 work quite well 
(Ref. 4], but the most commonly uséd and the one that is 
consistent with the potential enerjyy minimization formula~ 
fon 2S the strain snergy density [Ref. 2}. Because its 
employment reguires only minor software changes and it has 
been found to produce excellent results, #his refinement 
criterion was used extensively throughout the course of this 


esearch. 


Be. GRID ITERATION 


Another, father basic but less direc:é, method of 
Mecating regions of the mesh to bs refined is known as grid 
lteration, which can be implemented in one of two ways. An 


hitial coarse grid analysis may b2 repeated using Gither 4a 


j-4- 


1) 


Finer or a higher order nesh, and a2 comparison cf the resul- 


tants of interest between the two solutions will identify 
those areas of the demain where cafinement is most effec- 
tive. Another approach is based on the assumption that the 
greatest benefit is to be gained by refinement in those 
regions where a small perturbation, like the introduction of 
a single additional degres-of-freedion, produces the gre2atest 
change in the solution or one of its resultant parameters. 
Since additional degrees-of-freedom would be expected to 


produce the greatest change in those regions where ch 


1 





desired tesponse varied most rapidly, refinements based on 
this method provide results very similar to those obtained 
using the sclution parameter variation criterion already 
discussed. The grid iteration method may at first appear to 
be more computationally expensive, but 1+ was developed to 
be most efficiently implemented employing a special family 
of gslements. These elements, called "hierarchical", possess 
some very desirable properties for this application and will 
be discussed in the next chapter. 


Cw. ELEMENT RESIDUALS 


The major drawback with refinement criteria based on 
solution parameter variations is that their applicability 
appears limited to elastostatic problems. For this reascn, 
several investigators have recently developed arid refine- 
ment criteria based on element residuals, which appear pren- 
ising for application to all types of finite element 
problems, including nonlinear analysis. The reason for this 
is that the residual has essentially the same meaning 
regardiess of the probiem type [ Ref. 5]. For eéxamnple, 
consider the governing differential equation, 

D jew. ] = 2 = 9 
defined cn scme domain, where Df j is a linear or nonlinear 
differential operator, and the dependent variable u and t 
non-homogeneous term f are both functions of the independe 


a 


n 
n 
variables. Let the {inite element approximation to the 
Solution of the differential equation be Ui 7 u. Appiying th 


‘D 


Mmeere=Sntial Operator to the approximation gives rise to ‘th 


m 


residual, which is defined as 
av) 
R= ote |) 


and is not idantically zero unless 


tty 


ct 
Ye 
(Db 
tt 
t? 
: 
rw) 
q4 
iD 
Mm 
t— 
(D 
3 
iD 
Le 
ct 
(nh 
O 
— 
fer 
| 


meen is exac*. 


ae 
re xy 





The key t> uSing the residual as a criterion for deter- 
Mining regions of the domain where? refinement is necessary 
is the local residual on the element level, which indicates 
the contribution of the element to the total error of “tha 
finite element approximation. Since the residual is a Ddoint- 
wise quantity, the useful measure of the element error 
contribution is a norm of the element residual, or the inte- 
gral over the element of the product of the residual and 
some weight functicn. Bie anesgracion gis most readily 
performed using numerical quadrature. The grid optimization 
strategy then becomes one of refining the mesh so as to 
equi-distribute the element residual norms, by forcing then 
to become smaller and more uniform over the domain. 

There are some drawbacks to the element r2sidual refine- 
ment criterion which have not yet been fully resolved, such 
as appropriate selections of th? residual norm and the 


Mengnis £unction, and in the computation of she residua 
Soa 


2 a 


dtself. While the evaluation of th? rasidual oresen 


t- 
ct 


particular difficulties in the interior of the element, 


i 


is rarely determinable at the edges. The reason for this i 
that in formulating the finite alement approximation the 
element shape functions are defined so as to provide only 
that degree of continuity i to adequately model the 
physical problem; the most frequent consequence being that 
D Cu] is undefined along the interelement boundaries. 
Merortunately, this singularity canaot be ignored ard a more 
complicated analysis must be applied in order to bound the 


Pesidual contributions at these boundaries [Ref. 6]. 


De. A~POSTERIORI ERROR ESTIMATES 


A Sophisticaced extension of the element residual 
criterion is one based on computable error estimates from an 
initial finite element analysis. Thas utilizes che eneraqy 


a 





MOrmeGrncre Lesidual, in wWwoilech case =he weight function is 
the residual itself. Research in reliable error estimates 
was pioneered by Babuska [ Ref. 7, 8] for linear quadrilat- 
eral elements and more recently by Zienkiewicz { Ref. 9] for 
higher order elements. The major difference from the 
residual criterion is that instead of equi-distributing the 
element residual norms over the domain, they are normalized 
to compute error indicators for the elements, which are in 
turn used to compute reliable pointwise error estimates for 
the solution as well as the energy error over the domain. 
These quantities are of ovrimary importance because they 
provide not only an indication of where additional refine- 
ment is most effective, but also a measure or the quality of 
the mesh to determine whether additional refinemen* is 
necessary [Ref. 9]. Maemo pealewati ol strategy is to obtain 
meme arly Uniform, distribution of the error indicators 
throughout the domain, which corresponds to minimizing the 
error in the energy norm. The refinement procedure may prog- 
ress until all the lcecal errors are within some prespecified 
Pomerance. White the practical utility of such a refinement 
criterion is obvious, the mathematical devslopment and the 
algorithms involved are zather complicaced. However, tne 
process is not computationally expensive, and there now 
exists a prototype self-adaptive finite element code, FEARS, 


which implements this refinement criterion [Ref. 10]. 
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III. METHODS OF GRID REFINEMENT 
Once it has been determined where the initial element 
grid is inadequate and needs refinement, the next co 
tion is how the idealization in these areas shou 
improved. The choice of the refinement method to m 
may well be a more important decision than the selecti 
one of the refinement criteria previously discussed, since 
az least cne investigator has observed that for a particular 
method of grii refinement, the various refinement crit 
produce essentially the same solution resuits [Ref. 11}. 
ead beranemdient LS hs process of introducing edd2tional 
Gdegrees-of-freedom into selected regions of the finite 


element grid, and may be performed by one of three methods: 


(1) The pelynomial deaqree of the elements remains fixed, 
usually ata low value, while the size of th 
elements is reduced. This has become known as th 
h-version since element size is commonly denoted b 


mye £etter h-« 


(2) The size cf the elements, usually few in number, 
remain fixed while «he polynomial degree of the 
elements is increased. This has become known as the 
p-version Since polynomial ordjier is commonly denoted 


by the letter p. 


(3) The size of the slements may be reduced concurrently 
Wetman increase in theiz polynomial order. This is 
KnNOWn aS the cecmbined h- and p-version of the = 


element methed. 


is, 





A. CONVERGENCE OF GRID REFINEMENT 


It is well known that tha finicts e2lement xneathod 
converges with an increasing number of degrees-o 
an fact, fas as Che Gistltrcation for it de 
Therefore, the appropriate measure of the effectiv 
particular grid refinement method should be the associated 
rate ef convergence, which generally will be affected by the 
smoothness of the approximated function over the subdomain 
of interest. It has been demonstrated that wnen the refine- 
ment is performed uniformly over the domain, the associated 
Tate of convergence is asymptotic, provided the number of 
degrees-of-freedom is sufficiently large [Ref. 1]. The 
aSymptotic rate 9f convergence is often méasured as the 
Slope of the error versus cost curve in the linear, os 
Beymptotic, range when plotted on a log-log scale. Hg 
performing such an error analysis for the d 
Boemalaticn of the finite element msthod, EBhe 2E2OL 
usually the relative strain energy erro 
the energy norm, and the cost is assumed to 
function of the averages element size 
degrees-cf-freedom [{Ref. 12: p. 726]. Only in ¢ 
several years has there been any significant resear 
comparing the relative merits of the differant methods of 
grid refinemenc. Samecemen= SsOluUttons Of Elliptec boundary 
value problems are usually very smooth over convex domains 
except in the vicinity of corners, most of this research has 
focused on soluticns exhibiting SlIngularic.es, which 
severely hinder the rate of convergence, as in prebiens of 
MeeceUre mechanics and in domains with sce-sntrant corners 
feet. 1, 13, 14, 15). 

Mimccde=s Er a Einite element analysis to be both effi- 
cient and reliable, the asymptotic convergence zange should 


Fe entered for as few degre2es-of-freedom as reasonadly 
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possible. In general the p-version satisfies this require- 
ment better than the h-version. While it has been estab- 
lished that p-cecnvergence will nacessarily occur whenever 
h-convergence occurs, the converse is not true. For example, 
if the h-version using a uniform grid of linear elements is 
applied to a nearly degenerate problem, «he number of 
degrees-of-freedcm required for entry into the asymptotic 
range may be beyond the computer's round-off limitations, in 
which case convergence will not ccsur unless the polynomial 
ord¢r is increased [Ref. 1]. Numecical experiments on such 


a 


problems clearly indicate that +he p-version requires 
considerably fewer degrees-of-freedom than the h-versicn to 
achieve the same degree of accuracy. Recent analyses 
maets 1, $3] 5£ the asymptotic rats of convergence in energy 
for non-smooth solutions, “Using Uniform sefinenenc shia) 9 
sufficiently high numbers of degress-of-freedom, have demon- 
Strated that the p-version cannot have a slower rate of 


convergence than the a-version. Furthermore, if ¢ 


7) 
> 
Wy 


iee2tyY is confined to element boundaries, a 
case, «he error for p-mathod is invarsely prop 
the number of degrees-or- freedon, whereas the error is 
mayecrsely proportional to the squar2a root of the number of 
degrees-of-fraedom in the h-version. In other words, for 
this special class of problen, “H=meatS Os eon VELrqCenc] for 
the p-version 1S twice that for th2 h-version, which is dué 
Permarily to the ability of higher ordér polynomials to 
NMabsorb" singularities occurring at the element boundariés. 


Mes implies, at ieast for this typ2 of problem, <that i: 


Sy 


I-4) 


Seager ctO Minimize the efror for specified number 


a fo) 
degrees-of-freedom, the best str oe Noeno.) <0 subdivid 
ast o) 


e 
element of 


if 


2 
mae domain uniformly, but to use instead a siagl 
4 


increasing polynomial order [Ref. 15]. 
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cmcemite ts UNnLEKely that Ohn2 would attempt to solve 
such a problem using uniformly finsr grids, a more useful 
comparison between the convergence rates of the two versions 
would be based on adaptive refinement employing one of the 
solution-based criteria discussed in the previous chapter. 
It so harpens that the h-version, when used with optimally 
refined meshes, can have a higher convergence rate than the 
uniformly distributed p-version, provided that the element 
erder is sufficiently high. However, the p-version can also 
be employed with an optimal refinement criterion. While 
there are vet no proven theorems concerning the convergence 
rates for non-uniform refinement, Optaining the desired 
Semiyclon accuracy with Optimal p-distributions appears to be 
much less sensitive te the particular choice of the elements 
to be refined than with optimal h-refinement [Ref. 13]. 

It would ssem plausible that an even better optimization 
strategy would involve a proper combination of both the h- 
and p-versicns. It has been demonstrated for problems wit 
corner Singularities, ehhaw a prover sequence ar 
h-refinements combined concurrently with the proper sequence 
6 p-distributions results in extremely high convergence 
Mates, conjectured tc be exponential [Ref. 15]. However, 
moms DLOper combinaticn is difficult to determine, and adap- 
tive refinement based on the combined h- and p-versions 
poses some difficult data management problems. To avoid this 
problem a more promising approach, proposed by Babuska and 
Szabo [Ref. 1], empiovys @ graded nesh in which the element 
Sizes are first reduced according t>5 a geometric progression 
towards the singularity, followed py determining the optinal 
Beeesctcoibution for those elements using an adaptive 
Criterion. However, obtaining the optimal combination when 
employing this scheme can be &@ delicate matter and, astcund- 
ingly, the highest accuracy is achisved when the volyneomial 
Order of the elements actually increases with distance fron 


mie Singularity. 


a 
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There are some additional advantages of the o-version 
worth mentioning. Because the p-version employs fewer 


elements, there are lesser preprocessina and postprocessing 


Lf 


costs than for the h-version. Furthermcre, when bandwidth 
minimization and sparse matrix solution techniques are used, 
the solution time for the p-version is approximately the 
same as for the h-version for a specified number of 
degrees-of-freedom, and the p-version appears less suscep- 
Peepie to round-off errors. Finally, the p-version is simpler 
to implement adaptively than the h-version when hierarchical 
elements are 2mployed [Ref. 13]. 


Be. HIERARCHICAL FINITE ELEMENTS 


The hierarchical concept was first introduced as a 
Simple method for implementing the p-vérsion and as 32 
convenient device ‘fcr imposing boundary continuity between 
elements of different polynomial order { Ref. 9]. Since «hen 
a useful family of clements based on the hierarchical 
concert has been developed and incorporated into COMET-X, an 
experimental finite element ccde, disveloped by Szabo, which 
self-adaptively emplcys both the h- and p-versicns of the 
finite element methed [Ref. 14]. 

Homma =.C: GescriDrien of the hierarchical concept 
@onsider the conventional fin: element formulaticn which 


BPeoauces the follewing system of eguations: 


where mn is the number cf degrees-of-freedon, X(n) LS Ste 

Pex 2 global stiffness matrix, a(n) as fie, ELNS Et 

memsOoximaticn of the exact solution, and £69) 4 
O 


NM-component giobal icad vector. When n RE GRe:s VOEdce 
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degrees-of-fzeedom are added to the original systen using 
conventicnal refinement methods the system of Equations 


becomes: 


aa) (n+m) a 
Kn) a £ (Pare oe2) 


where the contributions to Ki...) aad RH each the refined 
elements result in a completely different set of coeffi- 
cients. If, on the other hand, *his refinement had been made 


hierarchically, the equations would become: 


K K 2 (M) sve 
a(n) a, (n Tm) a, av, 
= (Ci Cine. 5) 
K K ym) <(m) 
(mM, n) ~y(™m,m) 2° Y 


where K (n) and fy) are the stifiness matrix and force vector 
from the criginai system of equations for nh degqrees-of- 
Becedom appearing in Equaticn 3.1. However, Mis RO = the 
modal values of the finite element sclution for the mn 
additional degre¢s-of-freedon, but instead represents the 
Gifference between those values and the pointwise values 
obtained from the lower order polynomial interpolation for 
the unrefined mesh of n degrees-of-freedon. 

The primary advantage 9f hierarchical elements is imme- 
diately observable from Equation 3.3. Because the shape 
functions of an element of order Dp consestute a2 subsec of 
the shape functicns of an element of order ptl, the local 
PeeemnesSS Matrix and force vector for each hierarchical 
element is embedded in the stiffness matrices and force 
Beerors of all hierarchical elements cf higher order. 
Mrerefore, the global stiffness matrix ‘aly Se force vector 


ig ae | . , 
a eof the original system area preserved, thus saving 
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Sonscederaple Eim2 and E€Lfort expend2d on computing the coef- 
ficients fcr successive refinements [{Ref. 14]. Ancther 
advantage is that the hierarchical form of the global s*iff- 
ness matrix is more diagonally dominant than the one 


resulting from a conventional refinemert, resulting in 
improved conditioning and faster convergence when iterative 
solvers are employed [Ref. 9]. Another benefit of hierar- 
Chical elements, which arises from the “add-on" nature of 


the nodal variabies cf the higher order degrees-of-freedcn, 
Bempeia=) tne problemrof Maintaining boundary continuity 
between elements of different polynomial order becomes 
Peevial. Erstead of intrcducing global constraint equations 
for the higher order degrees-cf-fresdom, the nodal variables 
are simply set equal to zero and cond2nsed out, as if they 
were zero-valued Dirichlet boundary conditions [Ref. 2]. 
There are two major drawbacks with hierarchical elements 
that have hampered their widespread acceptability. The 
Sst, which has already been mentioned, 1s that the nodal 
Variables for the hicher order degrees-cf-freedom represent 
difference values rather than the nore easily identifiable 
values of the dependent variable itself. Secondly, when 
implementing the h-version cf tha finite element mnethcd, 
special integration rules must be introduced when the subdi- 
vided element 1s in hierarchical form [Ref. 9}. Q£ Course, 
the latter problem can be evaded by using the p-version, fer 
which the hierarchical concept was developed. In spite of 
the disadvantages of hierarchical elements, their consider- 
Person cOMputational efficiency and utility for grid optimiza- 
Seon wall certainly result in theic widespread utilization 


in future adaptive finite element sodés. 





Onee the analyst has identified where the initial grid 
needs enrichment and decided which refinement method to 
employ, he must then determine a systematic perecedure, or 
eegorithn, to perform the refinement according +t0 the 
criterion selected. The ultimate goal of such 2 prececure is 
t¢ design an element grid which neets the optimization 
objective of obtaining maximum solution accuracy fo 
ified analysis cost. While the analyst may or may not have 
pumerndicaticn of the accuracy of the solutior, he should 
Mave a Dieconceived notion of cost, or how much ertfert ke is 
Malling to expend to arrive at a better solution. Therefore, 
With some knowledge of the qdrid optimization techniques 
available and an understanding of the advantéeges ane 
vantages of each, the analyst can realize che grid optimiza- 
tion goal. 

There are essentially two adaptive grid optimization 


strategies: 


(1) Grid refinement, ween t hs initial  anelysis is 
performed on a relatively coarse grid, and new 
degre¢es-of-freedom are added t9 the same grid by the 
iterative ap Cilica.2.0n of the soluztion-based 


Grerrerzon< 


Wye Grid moditicaticn, in which the initial analysis is 
D 


th 


erftormed using a prespecified number } 
degrees-of-freedom, and the solution-based critericn 
ls employed to shift degrees-of-tEreedom from certain 
regions *o others. This may involve complete grid 
Eedecfinition in an effort to obtain a near-optinun 


meds an a Single cycle. 
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Much of the interest lately has been in the develooment 
cf cemplicated self-adaptive software packages which miri- 
mize the impact of the user's skill on the final solution. 
Ideally, the analyst would merely define the problem and ths 
program would automatically generate and analyze the optinun 
grid employing one or more of thes2 techniques, possibly at 


the user's cption. 


A. MATHEMATICAL PROGRAMMING 


No discussion of grid optimization techniques would be 
complete without a brief description of mathematical 
BaegqGauming, not only because it is how grid optimization 
was earliest attempted, Dias ere LMpOLcantiy, anes aes 
precisely what the engineer envisions when he hears the term 
Saacagazataon". It 15 not a grid optimization technique, fer 
se, but rather a numerical process of achieving any optini- 
zation obtective, once it is explicitly defined in mathemat- 
ical terms. In solid mechanics the finite element method is 
a numerical method for minimizing the potential energy func- 


tional, which in discretized form nay be written: 
Ti = 's u! kK =) eee Pane 4.1) 
= GG eR: ay 


where: 1s the global displacement vector 
+s the glicbal stiffness matrix, and 


2s the global lead vector 


In the classical finite element formulation, ‘the potential 
energy is minimized with respect to tha nodal displacements, 
meeen implies satisfaction of the following stationary 


Bondi * ions: 


of 





s71 
Daa = 0 (2 = 1,2,...,0) (Sane 4.2) 


where n is th? number of degrees-or-freedom. This leads to 
the very familiar system of linear equations: 


= 0 (EG. 4.3) 


CA 
C& 
i 

Crh 


However, since K and i are functions of the nodal coordi- 
Maces, then the nace energy could be minimized wit 
respect to the nodal coordinates as well. This would requires 
satisfaction Oe the following additional Stationary 


monditions: 
Ix, = O (ile eee 0) (Eqn. 4.4) 


where m is the number of nodal coordinates, Hac This aiffer- 
2 ae 


entiation leads to the less familiar system of non-linear 
equations: 
aK je! 
yg A a, (Gale? nsacyn) (Egn. 4.5) 
ay 9X4 Ox % 


mils, then, is the mathematical statement of the grid epti- 
Mization problem for the elastostatic cas The nodal 

displacement variables may be eliminated by pre the 
potential energy with respect to the nodal coordinates only, 
emogect +O the implicit constraint that Equation 4.3 is 
always satisfied [Ref. 4]. Unfortunatsiy this does not heip 
mich because the objective function is still nonlinear, 
rendering mest numerical optimization aigorithms inefficient 


and unreliable. The difficultv is even furcher compounded by 
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the requirement that the nodal variables be subject to sids 
constraints in order to maintain the defined boundary of the 
domain and+t> ensure that the elements neither distort 
excessively nor overlap one ancth2r. For all except the 
simplest of problems, ‘these constraints may be even nore 
severely nonlinear than the objective function, resulting in 
the analysis becoming prohibitively expensive [Ref. 2]. For 
this reason, nathematical programming in finite element grid 
optimization has been abandoned in favor of some equaily 
reliable, yet far mere computationally efficient grid opti- 
mization techniques. However, these early efforts with 
Mathematical programming were not totally in vain because 


they gave rise to the contouring techniques. 


Be. CONTOURING 


Since mathematical programming is infeasible for grid 
optimization, further investigations were conducted to 


suggest some guidelines t9 enable the analyst to construct 2 


gold wit Similar topological features of the numerically 
Smee m2 zed gric Without the computational ertort involved. 
maecke [Ref. 4], in employing mathematical programming in 


the soluticn of some sinple two-dimensional elastcstatic 
problems, observed that there was a very definite element 
pattern common among problems involving hiqh strain gradi- 
ents and that the nedes of the numerically optimized grid 
generally tended to be aligned along contours of some 
response function being modeled. Consequently, in periorming 
Smatyses on grids ‘whose construction was based on contours 
derived from an initial analysis, it was determined that the 
M@eiowing provided grid characteristics in regions of high 
strain gradients similar to the numerically optimized grid, 


= 


Pew as G4 Claction of the computational expense 


Pass, 





e contcurs of displacement 

e contours 9£ maximum principal stress 

e contours 9f maximum shear stress 

e contours of strain energy density (isoenergetics) 


e principal stress trajectories (isostatics) 


Since the strain energy density is the response which is 
consistent with the principle of minimum potential energy, 
isoenergéetics are the most commonly used contours along 
which element edges are aligned [Raf. 4}. However, there 
still remains the question of how to position the nodes 
along the contours. Fer this reason, isostatics have beccme 
increasingly popular because the ovrincipal stress trajecto- 
ties form a "flow net" of orthogonal curves which can guide 
the analyst in the layout of the elements [Ref. 16]. 

Saice GON°OULL 2G involves the redefinition of the grid, 
as opposed to a grid refinement, its primary advantage is 


in 
miteae che enriched mesh is not constrained +9 the elenent 


O 


configuraticn of the previous mesh. Therefore, there is a 
femee tO the amount of enrichment per cycle wnich can be 
performed and it is conceivable that an optinun mesh could 
be generated in a single cycle {Ref. 2]. However, while the 
computational cost of repeated analyses is reduced, the 
mee Orocessing Costs involved in constructing the contours 
and redefining the mesh c2n become quite high, ¢specialiy if 


the contours are complex. Aigorithms for performing these 


tasks in two-dimensional domains have been proposed 
Maer. 4&4, 17], but they are not extendabie to <~hree- 
dimensional problems. The major obstacle for two- and 


Maeee-dimensional doflains is that it is often difficult +o 
constrain the element edges to the contours without the 
elements becoming elengated or distorted to the degree that 
Mumerical inaccuracies result. Another difficulty, not 


ecdressed in the literature, is how the contour increments 
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should be selected when the response function is non- 


monotonic over the demain. 


C. SELECTIVE REFINEMENT 


The mest commonly employed grid optimization zechnique 
is that cf selective refinement. As its name implies, 
selected elements from a given mesh are enriched while the 
original element grid remains e2ssentially intact. The 
elements selected for refinement are determined by the iter- 
Meave aptlication of the solution-based criterion to indi- 
cate which elements contribute most to the solution error. 
The refinement can be performed by either the h-version or 
the p-version, or even the combined version if so desired, 
but the choic® is most often predetermined py the cavabili- 
ties of the available preprecessor. Since the addition of 
new degrees-of-freedeom over several iterations can guickly 
enlarge the problem, it is advisable to perform the initial 
analysis with a reasonably coarse grid of optimaily shaped 
elements, that is nearly square quadrilaterals or néarly 
equilateral triangles. Titcomis Cspecially important in che 
h-version where it is desirable +5 prevent the successive 
subdivision o2f elements from producing elongated new 
elements. One refinement technigue which will ensure this is 
mie sO Called “#ather-to-four scans" subdivision scheme in 
which a single quadrilateral or triangular elemenc is 


replaced by four new ones by adding and connecting midsid= 


nodes onthe edges cf the original eléemen= as shewn in 
mogure 4.1. The major difficulty in selective refinement 
arises when the addition of a node along an edge of the 


element to be subdivided creates a higher polynomial ordered 
edge for an adjacent element which is not te be subdivided. 
There results an incompatibility in the interpolation of the 


dependent vatiabie along this interelement boundary. Such is 
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Figure 4.1 Seme h-Version Subdivision Schemes. 


the case in the h-version scheme of Figure 4.1 and it also 
arises in the p-versicn when two ¢clements cf differe 

nomial order share a common edg2. When czhis situ 
occurs, the additional degrees-of-freedom do not actually 
represent degrees-of-fireedom act all because they must be 
humerically constrained to the polynomial interpolant of the 
lower order. Such ccnstraints are usually imposed in one of 
Basee Ways: global censtralnt equations ma 
constraints may be incorporated in the alen 
hierarchical forms may be used with the excess degrees-of- 
Ret. 2]. 


freedcm simply set to zero and condensed our [ 
nement techniques 


There are some cther selective refinen 
Meech GO Mot feqgquire any major software modifications. In 
the h-version, the continuity problem may be circumvented by 
employing any of the coarse-to-fine mesh transition schemes 
for which all of the element edges remain of the same poly- 
Bemeai order [Ref. 17: p. 210]. 3h2 


W 
to emelcy these schemes without permitting some elamenr 
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fas OLLLON, and the refinement must nearly always be 
performed interactively tather than automatically. For the 
p-version, interelement continuity can be easily ensured by 
employing varciable-noded isoparametric elements, Wee fi 
permit a single element to possess edges of differing poly- 
nomial orders [Ref. 7: p. 125}. 

The analyst must also exercise care when adiing new 
nodes tothe boundary of the domain to ensure that the 
appropriate boundary conditions are determined and applied. 
Micthermore, i2£ the bcundary is curved, the coordinates of 
the new node should be computed sSucno that it is placed on 
the actual boundary and not necessarily on the e¢4g 
element being refined [Ref. 2]. 

The important advantage of the selective refine 
mechnigue is that once an appropriate refinement cr 
has been determined, selecting candidate eslements for 
refinement in s¢ach cycle becomes straightforward. The 
refinement can then be continued indefinitely to achieve 
femiy high accuracy, but becaus2 the solucion phase is 
repeated for 2ach cycle, it is desirable to hold thr¢ number 
of cycles tc 2 minimum. Because th? nodes from the previcus 
mesh remain fixed for each cycle, selective “rs 


f n 

ideally suited for iterative solution methods. The so 

values obtained from the previous cycle, comnb qd 
m 


interpolated values for tha new degraes-of-fraeedom, orovide 
eameexcellent initial guess for the next cycle (Ref. 2 

The major disadvantage is that the limited amo 
refinement which can be verformed in each cycle may nan 
Meee several sycles te obtain an optimum grid. In addition, 
iz new degrees-of-freedom require interelement continuity 
constraints, data Management can become cumbersome unless 


the constraint is performed hierarchically. 
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D. SUBDOMAIN ISOLATION 


One of the obvious disadvantages of the selective 
refinement technigue is that the solution must be compietely 
mepeated £EOr 6€ach cycle when, in fact, fhe Tunber. os 
degrees-cf-freedom added in each cycle nay be few in compar- 
meanmce tne tocal for the problem. In addition, the number 
of elements requiring refinement in each cycle may only 
account for a small portion of the domain. Although the 
refinement criterion has indicated whare the grid is inade- 
Bate and the approximation is likely =o be poor, the solu- 
tion is repeated in each cycle for those nodes where the 
error is presumably small. Besides the apparent compu*a- 
mB=onat inefficiencies, this shortcoming severely restricts 
the amcunt of refinement which can be performed in the 
mIeEreglons of interest since it is desirable tS confine the 
size of the problem within reasonabie limits. An alternative 
approach is to reformulate the problem for those subregions 
where refinement is necessary and to accept the results of 
the initial analysis as an adequate solution fer the 
remainder of the domain. The elements requiring further 
refinement, Which ccnstitute isolated subdomains of the 
Se2ginal preblen, can generally be subjected to signifi- 
cantly greater refinement than would otherwise be practical. 
The solution obtained from the initial analysis is then used 
in imposing boundary valines on those degrees-of-freedom 
lecated on the boundaries of an individual subdcmain. These 
can, in turn, be used to generate the boundary conditions 
for any additional boundary degrees-of-freedom introduced by 


the refinement using an appropriate interpolation scheme. 


[res Gold SDtimization technique, which «he author terns 
"Subdcmain isolation", has some further advantages over 
selective refinement. The subdomain nay be selected arbi- 


Meaeaiy smati such that excellent results may be obtained 
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imoeaesinole Cycite using Wnhitorm refinement. Therefore, the 

Gesttveulcoaes sanvolvying cOoabpse-to-fins transition schemes, 
element elongation and interelement continuity can be 
avoided. Furthermore, one can choos2 as many subregions for 
refinement as desired without creating an excessively large 
problem. 

The obstacle which may prevent this technigue from being 
readily accepted is the notion that, by imposing erroneous 
boundary conditions on the subdomains, the convergence of 
the finite element method to the #2xact solution in these 
regions has somehow been tampered with. This aversion nay be 
somewhat abated by considering 4 Simple extension of 
Samnt-Venant'’s Principle { Ref. 18: p. 33]. Although the 
Sevaitions are not tigorsusly satisfied at the boundary, 
which may result in significant changes inthe response 
locally, the effect at some sufficient distance away will be 
negligible. The numerical evidence supports this premises. 
While errers in the boundary valu2s may soméwhart restricc 
the accuracy 9f the dependent variable, great improvements 
can be realized in the accuracy of its gradiants, which is 
Gere Otten the goal cf the optimization. Since the initial 
analysis provides the boundary valuss for the subdomains, ict 
is desirable that 1ts solution be as accurate as teasonably 
possible. Fortunately, since subsequent refinements are not 
Perrtormed on the original grid, che initial analysis may 
involve significantly mor2 degrees-of-freedom than in the 
case of selective refinement. 


Ee MESH GRADING 


Greeetial Grid Optimization technique +o be discussed 
employs a mesh for which the element sizes are successively 
reduced, according +o some geometric sequence, towards a 
Selected region of the domain. Gio taaiee erg that nesh 
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Taadimgeis mot really an Sptimization technique sinc? it is 


most often applied on an "‘a-priori"™ basis rather than acap- 
tively, and that it does not lend itself well to the itera- 
tive application of a solution-based criterion. However, the 


er 
technique is Simple to use and its implementaticn requires 
few software modifications. Furthermors, a solution-based 
refinement criterion can be used t> give a measure of the 
quality of the mesh to indicate wasther a more pronounced 
grading may prove beneficial. Depending on the solution 
parameter of interest, mesh grading can provide excellent 
accuracy at a low analysis cost. This refinement method must 
therefore be considered among the grid Spe a Zee Lon 
techniques. 

For the less elaborate finite element preprocessors, 
mesh grading is often the only refinement means avaiiable 
without rescrting to a uniformly finer mesh involving many 
more degrees-of-freedcm. The most common méthod of implemen- 
tation in two-dimensions is to first d¢fine the probilen 
domain in terms of a curvilinear quadrilateral by seiecting 
four keynodes along the problem boundary. Then the boundary 
nodes are spaced according to som? geometric sequence based 
on the user-provided bias parameters which determine the 
density or the nodes towards selacted points on the four 
quadrilateral edges. Finally, curves are generated +o 
connect tke boundary nodes on oppesits edges of the gquadri- 
lateral, thus producing a graded mash. This process, which 
Can also be extended to three-dimnensions, is the nesh géner- 
ation scheme employed in the finite element cede GIFTS 
[Ref. 19]. 

The major disadvantage of mesh grading is that in order 
to achieve sufficiently small elements in tha regic 
interest, the e2lements must grow successively larger away 
Bom that region. This May be very undesirable, eésoecially 


if refinement is called for in mors chan one region of the 
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domain, in which case the mesh must be generated and gradea 
by subdomains, thereby complicating the data management 
involved. 

Another disadvantage is that unless the demain possesses 


some special geometric symmetry, axcessive element 


Figure 4.2 Graded Mesh for a Perforated Square Plate. 


elongation will usually result if a highly prenounced 
Meadang 1S required. Some configurations are particularly 
well suited for refinement for mesh grading such a the 


W 


S 
classicai perforated square plate oroblem in solid nechanic 


shown in Figure 4.2. 
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Now that the necessary tools for performing grid cptzimi- 
zation have been intrecduced, it is time to employ them in an 
attempt +o obtain optimal solutions to some practical prob- 
lems in engineering. An obvious starting point for such an 
investigation is the one-dimensional boundary value probien. 
Watle most cf the fruitful research in grid optimization has 
concentrated on problems of higher dimensions, the one- 
diménsional problem is a very convenient device for studying 
Finite element grid optimization. Foremost, one-dimensional 
finite element models possess a unigue connectivity in that 
adjacent elements meet at their end nodal points. Therefore, 
refinement by the h- or p-versions, or by relocating nodal 
points béccmes a trivial task, which does not involve any of 
the difficuities so frequently encountered witn higher 
dimensicnal problems, such as preserving interelement conti- 
nui+y and maintaining optimal element shape Furthermore, 
one-dimensional studies can often provid? valuabie insight 
Bom the solution of Mere Girticult higher dimensional 
problems. 

The primary concerns in the selection of the preblems to 


Peestuaqied were as fcilows: 


(1) there should exist an analytical solution to provide 


@améeans of reliable error analysis; 


Mee cne solution and its resultants should exhibit 
sufficiently high gradients so that the effective- 


ness cf the grid cptimization is readily observed. 
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several solution paramet2rs which are easily con ; 


ute 
requiring minimal software changes to an existing finite 
element code. 

Furthermore, for the one-dimensional investigation, it 
was decided to simplify the analysis by exploiting the 
linear sglements. While it is grant2d that improved solution 
accuracy may generally be obtained by employing higher order 
elements, it will be assumed that conclusions based on the 
use of linear elements can be applied as well to elements of 


higher pcelynomial order. 


Ae. ELASTIC CABLE PROBLEM 


Consider an elastic cable under tension Tf, stzetched 
between two points a distance 2L apart. WE ene. Cable ws 
sum@ported DY 2 Winkler, or elastic, foundation of modulus k, 
and a concentrated load P is applied at the midspan, the 
resulting deflection v(x), (O75 x SS LL}, 25 as shown “in 


Beaure 5.1. Miew@earelyeteal Solution and finite el 


(D 


nent 
approach for this preblem are presented in Appendix A. 

The initial finite element analysis was performed using 
Peminear Clements of uniform length. From this initial 
analysis the approximate distribution of the following solv- 
tion resultants was cbtained over the domain: 


e the displacement, v (the solution) 
e the slo 


tn 
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e, v' 
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ct 
ty 


® che in energy, JU 
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ct 
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* the s in energy density, SED (dU/dx) 
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Figure 5.1 Yension Cable Deflection on Elastic Foundation. 


Subsequent analyses were performed for finite element nodeis 
using the same number of elements, but with the nodes redis- 
tributed to achieve approximately uniforn variation of the 
above parameters over each element. Note that the strain 
energy refinement criterion produc2s slements of identical 
Strain energy content. In addition, ‘the problem was solved 
employing graded element models of various adjacent element 
length ratios. The resulting elemeant models based on these 
refinement criteria are shown in Figure 5.2 (a-f). The 
graded model (b) fer an element langth ratio of 1.2 is 
presented for comparison because it prodtced good overall 
solution results. 

AS previously nentioned, wie esOMitaon PesuLtants oO 
primary interest are the maximum displacement, the naxinun 
Slope, and the integral of the displacement over the domain, 
because they represent important analogous soiution results 
in nearly all fields in which finite element anaiysis is 


orten performed. The accuracy obtained in these values for 
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Pigure 5.2 Tension Cable Refinemerts. 
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each of the refined models is presented in Table I. As ca 
Pemseenein ragvec 5.2,05 one Strain energy end strain energy 
density criteria preduced extrem variations of #?2i¢ 


m 
Jength while the criteria of displacement and slope res 


qn Gases “EEE ES ae ee me = a ee ee ae a ee 


: 


TABLE I 
Tension Cable Problem Solution Results 


Pe ep a | 


| 

| Prceblem Parameters: L = 100 in. 

K = 1 psi 

| T = 1000 Tp 
P = 1000 lb { 
| : 
Variation Percentage Relative Error 
i Refinement | 
| Gaoter on v(max) =v! (max) ix | 
! Uniform -0.40 0.36 2 | 
Graded (1. 2) -0.19 0.07 eas } 
v -0.18 0.06 0.39 | 
vi -0.23 205 Oeo7 | 
U -1.03 0. 05 Beis | 
S538) -1.29 02705 a.17 | 
| U (mo a) -0.53 0.04 1.63 
| 
i SED (mod) -0.51 03 1.48 | 
ee a ee 
in more moderate variations. It can be observed in Table I 


that the more pronounced refinements based on energy distri- 
PeerGn result in greater accuracy for the maximum siope but 
Maeh the accompanying severe penalty of significantly poor 
estimations of the maximum displacement and the ah eiee 
Stent icy. Lomi see lee CUlakr oprobiem the unitcrh ¢ 


provides optimal accuracy of the integral Giaa<eL ty, 





therefore refinement cannot reduc®= its error. Yet great 
improvement in the accuracy of ths maximum slope and modest 
improvement in the accuracy of the maximum displacement can 


be achieved with moderate refinements based on the dispiace- 


ment and slope distributions. 


’ 


One might assume, and correctiy so, that the ability of 
the energy refinements to produce the best accuracy for the 
Maximum slope is due to the extremely small elements which 
result in the area where that quantity occurs. Furthermore, 
it would be correct to propose tnat the reason for these 
refinements preducing poorer estimates than the uniform 
model for the other two quantities of interes* is that the 
excessively large elements in the regions of low yradiencs 
severely overstiffen the model there. Tt would then seen 
Peausazble tc improve the accuracy for the maximum displace- 
ment and the integrai quantity by redistributing «he nedes 
in these regions to rrevent such excessively large elements. 
This was done by arbitrarily employing a grading scheme to 
the three largest elements to produce the modified refine- 
ments based on strain energy and strain energy density shown 
in Figure 5.2 (g) and (h). AS can be seen in Table I, such a 
modification did indeed significantly reduce the errors in 
the maximum displacement and the integral, bux it even 
further improved the accuracy for the maximum sloce. 

One might conclude from Table I that thse best overall 
model was obtained using the graded mesh, and that since it 
is €asier to obtain, it should be deemed the optimal grid. 
Meee nts Particular grading was cnosen for precisely chat 
rsason and was presented only as a means cf comparison. [In 
practice, the selecticn of a grading ratio is somewhat arbi- 
trary and making an adequate choice nay be difficult. 


There is oe confusion as to which refinement 
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Peoduced the "best" solution scoeac? for this problem an 
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Pode Letinenent.  Betore any Optimization process can be 
pursued, the optimization goals must be explicitly defined. 
Clearly, as is the case in this problem, +he designation of 
the cptimum grid would depend haavily upon which of the 


three solution resultants is most critical to the analysis. 


Be. TAPERED BAR PROBLEM 


The linearly tapered bar under axial loading has 
received considerable attantion ani was cone of the early 
problems for Wieeheoanalytical Gia der Ot Liza 2 On was 
employed. Consider a tapered elastic bar of length L and 
modulus E, fixed at one end, with an axial load P applied at 
mie Other, . for which fhe axial displacement u(x), 
(Oo < x € Ll), is desired. The cross-sectional area varies 


linearly from A, at the tived en@ec] AU a= the tip, as snowa 


iG 
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Pigure 5.3 Linearly Tapered Bar Under Axial Loading. 


me Figure 5.3. The analytical solution and finite alenent 
approach are presented in Appendix 3B. 
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One "of them sigqmeficant features of the t 
problem is that the maximum stress can be very dif 


model accurat2ly, and it is for prsciseiv these problens 
exhibiting large strain gradients that grid optimization 
becomes most beneficial. Incvessstingly, +he stresses 


obtained at the element midpoints are exact for this 
problem, and the difficulty arises frem the inability of the 
constant slope shape functions to model the maximum stress 
occurring at the boundary. In examining this problen, 


Prager {Ref. 20] demonstrated analytically that when each 


ct 
oe 
ay 
iD 


element has the same strain energy content, the rela 
error in displacement is identical for all the nodés. 
However, this phenomenon appears peculiar to this preb 

and the author does not subscribe to such a méasure of an 
optimum grid. Judging =he effectiveness of a particular 
refinement based upon the deviation or the mean value of the 
pointwise errors generally tends to be unfavorable to cpti- 
mization procedures since the aimost always introduce many 
mor? nodes in those regions where the response is nest 
difficult te model. Hence, an improvad solution may have a 


larger mean value of the pointwise errors [Ref. 3]. 


The criteria employed in che refinement of the tapered 
bar med¢l are identical to chose used in the cable problem 
and their effects are shown in Figur2 5.4 (a-e). Two 2xcep- 
Meons are that nowth> displacement and strain ¢nergy 
criteria produce identical refinements, and the graded nodei 
chosen as the best overall is now basad on a cgcrading ratio 
@eni.4, producing a more drastic refinement than that of 1.2 
for the cable. tao pune ner accnmOnsckaces the difficulty 
eeyOlved in obtaining adequate clement grids ergl ely 


Ma prEiori basis. 
The Solution results are presanted in Table II and the 
MOSst readily apparent observation for this problem is the 


large errors in the maxinun slops, which would severely 
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Figure 5.4 Tapered Bar Refinements, 
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TABLE If 
Tapered Bar Problem Solution Results 


-_ } 
| 
| | 
| | 
| Problem Parameters: L = 100 in | 
| A = 10.5 in2 | 
i A = Ge > ine | 
B = 110x106 psi 
| P = 10x103 lb 
; | 
i Variation Percentage Relative Error | 
| Refinement L 
Gites .On u (max) ut (max) udXx 
| 
| Uniform =3250 =37.5 0.68 | 
| Graded (1.4) Silos ke: “4.1 0.1714 | 
ea -0.85 - 10.6 oP bs. | 
u? -1.81 -7.7 0233 
Sb -6.54 -3.6 Tere 
| ————— ee foo | 
SED (mod) -1.99 -3.6 O36 
| 
{ 


' 
| 
| 


underestimate the maximum stress. These results are based on 
quadratic extrapolaticn of the exact slopes at the element 
Medpoines, since the linear shap= functions would produce 
even poorer estimaticns of the maximum slope. As before, the 
more extreme refinement based on the strain energy density 
variation provides the most accurate estimation of the 
maximum slope, but with the accompanying degradation in 
estimates for maximum displacem2nt and the integral of 
displacement. Again, the large errors in these values my 
be Significantly reduced by ¢mploying a grading schame to 
restricz the size of the larger elements as shown in Figur: 
Seeeqat). Unlike the Efevious probléen, such a nodification 
has no effect on the estimate of naximum Slope because of 
the extrapolation of the slement midpoint slopes, which ars 


Smace S=acartdiess of the element modal. 
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A different version of the tapered bar probien, Rohm 
which the displacement and strain anergy criteria will not 
produce identical refinements, involves replacing the 
Bemcentrated tip load Pwith a linearly varying axially 
distributed load q(x), specified by the values at the fixed 
end q, and the tip q,. The problem may be further modified 
by reversing the bar such that the maximum slope occurs at 


the fixed end, while the maximum displacement occurs at the 


UX) 


ee ay 


Bogure 5.5 Reversed Tapered Bar with Distributed Load. 


free end as shown in Figure 5.5. The case of the linearly 
Merying distributed load is included in the formulation in 
Beeendix B. 

This problem was solved for a uniformly distributed load 
using tne same precedure 2s in the previous two problems. 
The refinement models are presented in Figure 5.6 and the 
solution results in Table III. The observations are consis- 
tent with those made in the previous problems, but now one 


would likely 2gree that the refinement based on the strain 
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Pigure 5.6 Reversed Tapered Bar Refinements. 
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TABLE {il 
Reversed Tapered Bar Solution Resuits 


ye ee 


| 
| 
| Problem Parameters: L = 100 in 
| A = Jie) she 
| A = 10.5 ine 
E = 10x106 psi. 

| q = 100 Lb/in 

Veaeoe con Percentage Relative Error 
Refinement 
| Chicos ron u(max) uf (max) gudx | 
| Une fo Sin SS. =O -7.1 | 
| u Bal) ares, =2 0 
| ut meee Ole =i | 
| U S25] ano “u.7 | 
i SED el} pee: -3.4 Sia 5-5 | 
| 
| SED (mod) = 310 -3.4 -4.0 
| 


energy density would represent an optimal grid, provided 


that modifications are introduced to prevent any elements 


from growing 2xcessively large. 


C. GUIDELINES FOR ONE-DIMENSIONAL GRID OPTIMIZATION 


The most important lesson to be learned from this one- 
dimensional study is that the grid optimization procedure is 
necessarily dictated by the optimization goal, or the under- 
lying purpose for performing the finite element analysis. No 
element grid can possibly provide optimum accuracy for every 
solution resultant cf interest. In solving these simple 
problems, a balance has been sought for achieving adequate 


eceuracy for three of the mor2 important solution 
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resultants, with emphasis on the maximum value of the deriv- 
ative of the dependent variable, Watch ne= «= Cheer LS net 
only the mcst important part of the solution but also the 
BoSt difficult to obtain accurately in finite element 
analysis. 

MicmerheOrtateGr2a Optimization techniques of intro- 
ducing mcre degrees-cf-freedom by subdividing the éelenents 
or increasing their polynomial order have been intentionally 
Si]atted in favor of the Optimization strategy of seeking 
Maximum solution accuracy for a specified number of 
degrees-of-freedom using linear ¢lements. This 1s because 
such a procedure is not so straignttorward iis wO- 
dimensional problems where the number of degrees-or-freedonm 
are dependent on some geometric considerations, which do not 
appear in problems of one-dimension. Based on this choice of 
optimization strategy, it appears the strain energy density 
variation provides the most useful criterion for the 
tive refinement of the initial grid. Yet all three pro 
demonstrated some pathological results that can arise when 
the elements are permitted to grow excessively large in the 
regions where the strain energy density varies the least. In 
applying a scheme tc restrict tha size of the largest 
elements, no nention has been made of how tec determine when 
an element is excessively large. Iz has becene the experi- 
ence of the author that any element representing over half 
of the demain should frobably be considered to0 large, and 
measures should be employed to restrict its size. 

Tt would appear, at least for these classes of probisnms, 
mm@ae this difficulty of decreasing accuracy of a particular 
solution parameter for successive refinements can be ignored 
by merely accepting the largest alu2 among the cycles as 
the mest accurate sclution result. For example, it was 


5 


demonstrated that the refinement based on strain energy 
density provided significant improvement in the accuracy for 


>| 





the maximum slope but underestimated the maximum displace- 
ment even mere than the initial uniform grid. Assuming that 
+he linear element mcecdel always underestimates such maxina, 
the maximum slope for refined grid and the maximum dispiace- 
ment for the unrefined grid could ba accepted as the optinal 
results of the analysis. The fallacies of such a procedure 
eee that, first, the refinement may not represent the 
optimal grid as it has been defined and, second, for self- 
adaptive finite element codes the user is provided with the 
feotinum grid" of the final cycle and the solution results 
memerecf. 

Melosh and Marcal [rRef. 21] have proposed an alternative 
use of the refinement criterion based on strain energy 
density variation which avoids the problem of excessive 
element rowth altcaether. Beginning with a reasonably 
coarse uniforna grid, those elements with the greatest strain 
energy density variation are selectivaly refined by either 
Siedavrding them Or increasing their polynomial order with 
the intreduction of additicnal dejgrees-of-freedon. While 
such a procedure does not equi-distribute? the element strain 
energy variations, it can reduce them all to some presseci- 
fied tolerance, such as 3 percentage of the average element 
variations for the initial analysis. Because this procedure 
is particularly attractive fer grid refinement in problems 
of higher diménsions, it will be employed extensively for 
the study of grid optimization for two-dimensional voroblems 


miecne next chapter. 
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VI. APPLICATION TO 


TWO-DIMENSIONAL PROBLEMS 

Since investigators began working in the field of init 
element grid optimization in the early 1970's, nearly allo 
the effort has been devoted to the development of a systen- 
atic procedure for obtaining optimal Geidss for *wo- 
dimensional problems of elasticity. Even todav there are 
several competing approaches to this problem and no partic- 
ular one has yet been overwhelningly accepted as the 
Peorerred method S£ grid sptimization. While it is tae two- 
dimensional problem for which most of these techniques have 
been developed, their application *9 such can be much more 
difficult than for the one-dimensional case. Almos* invari- 
ably when parforming grid refinement on two-dimensional 
domains, the analyst is confronted with the problems of 
Maintaining interelement compatibility and preventing severs 
element distortion. 

In selecting an appropriate two-dimensional problem for 
the application of seme grid optimization techniques anda 
comparison of their effectiveness, it is desirable that the 


test case possess the following properties: 


e the analytical solution should exist in order to 


perform reliable errer analysis 


Seche solution should exhibit SUtmucwentiy large 
Gradients to previd? a meaningful measure of the 


Betinement effectiveness 


= 


e the idealization should have cne degree-of-freedom per 
nodé and possess sinple boOUndasy —CONndi=Lens <0 
Minimize the computa tional arfort lnvoived ia 


repeated solutions 
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There are few problems that meat these criteria, Dut 
Banceavenanme ©OrSston Of & non-circular section oro 


good test case. 


A. PROBLEM DESCRIPTION 


Consider a solid circular shaft of radius "a" made fron 
isotropic material of shear modulus G and having a circular 


groove, cr keayway, of radius b along a generator of tke 


y 
° — xX 
Figure 6.1 Cross-section of Shaft with Keyway. 
Smee. Lhe shart cross-section 1S shown in Figure 6.1. The 


Shaft is subjected to an applied torque T which produces 
goangle of twist per unit length 9. The ad 
solved by finding the Prandtl torsional stre 


Ss u 
Which satisfies the governing differential equation: 
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BD CCuLON ET RemD=eernechlst condition that ¥Y = 0 on 
section boundary. The torsional stress function is dafine 
such that the shear stress T at any point on the domain ma 


be expressed 2S; 


tT = G0 [( B24 ak CSG ne O = 2) 


For this formulation, the angle of twist 9 is prescribed, 
Beaener than the applied torque fT. Rie coecler cmos Gal. cu 


lated from the area integral: 
T= 260] WV dA (Eqn. 6.3) 


mre. analytical soluticns of Equations 6.1 and 6.2 are 
derived by Sokoinikof£& [Ref. 22: pp. 141-143] ard are 
presented in Appendix cC along with the evaluation of 
Equation 6.3 anda prescribed finite element formulation. 
fee this problem, the three solution resultants of interest 


Boe the grid 2otcimization study are: 
(1) maximum value cf the dependent variabl¢e, or torsion 
Eunct aon Wee 


(2) maximum value cf the gradient of the dependent vari- 
able (a quantity proportional ‘*o maximum shearing 


stress Tmax); 


(3) the area integral of the dependent variable cver the 
domain (a quantity proportional to the applied 


EGE qWer). 


These quantities - the dependent variable, its gradient, and 
Ss 


an integral thereof - are selected as representative of 
entities whese error one might wish to ninimize in a finite 
element analysis. 





B. COMPUTER IMPLEMENTATION 


R= ceneoe seen 2a figure 6.1, the domain of this problen 
is symmetric about the x-axis, therefore the finite element 
solution need only be obtained for the upper half o£ the 
domain. For all of the solutions presented herein the 
problem geometry is defined by assigning the dimensionless 
ratio, b/fa = 90.4, and an acceptable upper limit on the anal- 
ysis ccest was arbitrarily chosen to be that corresponding to 
approximately one hundred nodal points. The computation and 
assembly of the finite element matrices and solution of the 
resulting system of equations was p3rformed using the steady 
Beaece heat Condwetion operations of CAL-NPS [Ref. 23]. This 
group of subroutines comprises an efficient finite element 
code fcr solving Poisson's equation in two or three dimen- 
Sions and has the additional advantage of overmitting 
variable-noded isoparametric elements. 

Since there was no readily available interactive prepr 


. 


cessor which lent itself well to adaptive aesh refinement, 
mye author had no choice but to create his own. Since the 
problem domain is simply connectsd, the automatic mesh 
generation was performed employing inverse epecnd) Cn wc 
Single cubic isoparametric element 9f the oe. =a 
Meo the problem dcmain [Ref. 24: pp. 228-229}. Ma 
boundary nodes were repositioned to conform to the actu 
domain beundary and additional nodes generated during 
refinement process were mapped using the same procedure. 
Since the finite element code salected for this investi- 
Secron provided output only for th2 nodal values of the 
dependent variable, it was coupled to the author's pos«pro- 
Gessor. Such a postprocessor is necessary in tke optimiza- 
tion process for computing nodal values of shear stresses 
and strain energy density, SVoteneconetnlputLons tO. Torque 


ogee -Ocal Strain energy, 2nd exact results from theory. 
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C. ASYMPTOTIC ERROR ANALYSIS FOR UNIFORA REFINEMENT 


The ccncept of asymptotic convergence rate for unitornaly 
refined grids was presented in Chapter 3. When the number of 
uniformly distributed degrees-of-freedom is sufficiently 
large, the log-log plot of the relative energy error versus 
the number of degrees-of-freedom is approximately linear in 
the asymptotic range. The slope of this line represents the 
asymptotic rate of convergence in 2nergy. 

It so happens that relative error in the torque T of 
mas prceblem is equal to the relative energy error and 
therefore exhibits this linsar asynuptotic behavior on the 
log-log plot against the number of uniformly distributed 
degrees-of-freedcm. FOSeiiacciV,ee cas Ovass two “SOLUTION 
resultants of interest behave similarly. This will prove 


very beneficial in performing the error analysis for this 
two-dimensional study for two reasons. First, because it is 
mimecessariliy difficult to construst an ot. Csi a eich 
the same number of degrees-of-freedom as a uniform grid, the 


linear behavior of the solution resultants in the asymptotic 
Banage On the log-log scale permits interpolaticn fcr any 
number of degqrees-of-freedom. Then the solution results for 
a uniform grid of the identical number of degrees-of-freedonm 
provides a reference fcr comparison to determine the effec- 
tiveness of the optimization technidque. Secondly, oii the 
Benvergence cate of a particular solution resultant is 
extremely slow, as is often the case for maximum stress, it 
Beeomes ditficult to gain an appreciation of the true e 
tiveness of the optimization. For example, an order of 
magnitude reduction in the relative solution error 
require an order of magnitud2 increase in the number of 
degrees-of-freedom using uniform refinement, but re 

few additicnal degrses-of-freedom using an opti 


Za 
Beenmique. Therefore, it will be enlightening t9 extrapolate 
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Figure 6.2 Uniitorm Linear and Quadratic Element Grids. 





Poeeitclative srrer versus degress-of-freedom curve te obtain 
a rough approximation of the number of degrees-of-freedom 


Mecessary to obtain solution accuracies similar ‘tc th 


iD 


th 


optimal grid, but using successively finer uniform grids. 0 
course, this is only an 2stimation and ignores such reali- 
ties as numerical ill-conditioning and computer rcund-off 
eaeror. 

AeuGLtOrm grad is e@ne fOr which all of the elements are 


of the same size h and the Same polynomial order p. Clearly, 


Meets impossible to cbhtain such a gri Ge chs pare ci let 
domain using isoparametric mapping, but a nsarly uniforn 
grid may be constructed in which the elements are of approx- 


imately the same size. Such uniforn grids ar? shown for the 
cases of linear quadrilateral elements and quadratic seren- 
dipity elements (Fig. 6.2). For this geometry, the uniforn 
grid is not uniguely defined for a specified number of 
elements. This is because, in performing isoparametzic 
Mapping, cher must be specified four keynodss on the actual 
@emain becundary to ccerrespond t¢6 the four corner nodes of 


the parent square. SinGe this “Gemean has only three 





Figure 6.3 Keynode Placement for Isoparametric Mapping. 


oP, 





vertices, the placement of the fourth keynode is a 
Mesctetion Ct she analyst (Fig. 6.3), and can have a notice- 
able effect on the solution results. 

The asymptotic error analysis was performed for the 
three solution resultants of interest uSing uniform grids of 


linear and quadratic elements. Th2 results are presented in 


( 


Figure 6.4. All of the solution resuitants behaved as 
predicted with the exception of the maximum torsion function 
value using linear elements, m(1). It appears that the 
Mecircracy of this particular parametsr is very strongly 
dependent upon the keyncde placement. The curve constructed 
in Figure 6.4 represents an average ‘fcr saveral keynode 
Bosicions. 

While Pigure 6.4 is intended primarily +9 serve asa 
reference tool for future analyses, it provides some inter- 


Pacing irformation: 


(7) For the cases of maximum ‘torsion function value and 
applied torque (and 32nergy), the asymptotic rate of 
convergence using quadratic alaments is more than 


*wice that for linear elements. 


(2) While the error in torque for the quadratic case is 
@lways smaller than that for tha linear cas2, the 
linear grid may provide better accuracy for the 
maximum torsion function value max in the opre- 


asymptotic range. 


(3) Both accuracy and convergenc? rate in the naxinun 
Shear stress are only minutely greater for the quad- 
Batic element grid than for the linear one. 


However, tor this last observation, the point must be nade 
that for the linear element grid, the mnaximum shear stress 
Was Obtained by quadratic interpolation rather «han from the 


iznear shape functions. While this will greatly imcrove the 
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accuracy cf the maximum shear stress approximation, it wiil 
have no effect on its rate of convergence. Pies =fon =, 24 
obtaining an »ptimal estimate of the maximum shear was the 
purpese of the analysis, there 1s much to be said on behalf 
of linear elements besides their computational efficiency. 


Of course, this observation is based on uniform grid refine- 


( 


ment, which would rarely compete favorably with the ovtimi- 
Zation technigues tc be examined. 
The reason that the rate of convergence in maximum shear 


stress is so poor usin uniform refinement tor this problen 


T CAD=0 


Figure 6.5 Stress Distribution on Shaft Keyway. 


can be seen in Figure 6.5. The sh22r stress varies greatly 
ieee ad ShOt. distance, by Lncreasimg from zero at point A to 
its maximum value at point B. As 3 Suc there exists a 


region of excessively large strain gradients 
keyway which severely hinders the rate of converge 
uniform grid is employed. If the keyway radius were allowed 


memaveroach Zero producing a Singularity in the solution, it 
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would likely be necessary to employ even higher crds 


Dp 


elements via the p-version in order to achieve convergenrc 


using uniform refinement (Ref. 1}. 


D. PROBLEM SOLUTION WITH GRID OPTINIZATION 


The finite element solution of the torsion problém will 
be obtained employing the following grid optimization tech- 


nigues as presented in Chapter 4: 


Seececntouring 
"MEOnNtOuts Of the tOrESiOn function: linear elements 
- contours of shear stress; linear elements 


- contours cf strain energy density; iinear elements 


e Selective Refinement 
- h-version; linear elements 
- h-version; quadratic elements 


- p-version 


eesubdcmain ILsolaticn 
- Jlinear elements 


- quadratic elements 


e Mesh Grading 
- linear elements 
- quadratic elements 


1. Contouring 


The original finite element analysis was performed 
cn a uniform grid of 98 linear elanents, 78 nodes, and 72 
degrees-of-freedom. The finite element soiution provided the 
Megat Values of the torsion function w , from which the 
conventionai nodal resultant ©6f shear stl2ss Tt ,; and 


applied tcrque fT were computed. Based upon the naxinum and 





Minimum values obtained for each parameter, alcng with 
Bamcolderat.on f96 thear Values along the boundary, The 
contours te be used for nodal placement in ach case were 
selected. The number of contours for each case was chosen *to 
Maintain approximately the same number of degrees-of-freedon 
as for the initial analysis. 

The points for each contour value selected were 
obtained by linearly interpolating between the nodal values 
of each parameter obtained from tha initial analysis. The 
contours were constructed by smoothly connecting the points 
by hand. The element layout along the contours posed the 
most formidable problem because ths coarse-to-fine tran- 
sition often resulted in severe elsment distortion, and it 
sometimes became necessary to degenerate quadriiateral 
elements into triangles when the transition was acut2. [It 
was decided chat the optimal element shapes should be 
preserved along she contours in regions of highest stress. 

The contours obtained and the corresponding grid are 


presented for each of the following solution resuitants: 


e torsion function (2g. 25.5.6) 
e shear stress CEs: D7} 


e strain energy density, SED (Fig. 6.8) 


The resulting grid for each of the response function 
contours produces smailer elements in the région of greatest 
stress near the bettom of the ke2yway and around the 
periphery of the shaft where the stress is moderately high. 
Consequently, the elements near the center where the stress 
is zero are larger. Mhesie, Of »course, are the desired 
effects for an optimization criterion. A somewhat unusual 
behavior is observed at the point of intersection of the 
keyway and the shart boundary wher? the stress is also zero. 


Apparently, the shear stress gradient is larger than the 


fu 


Sead2en=s in torsion function and the strain energy 


ensity, 
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resulting in smaller elements being produced in that region 
by the shear stress criterion (Fig. 6.7) than by the octher 
criteria. While all of the grids possess to Some extent tne 
desirable features cf an optimal grid, the strain energy 
density function produces a far more drastic refinement 
towards the point of maximum stress, while the others repre- 
sent more moderate refinements. In fact, the SED contours 
are sc dense around the keyway that the coarse-to-fine 
element transition scheme must include degenerate quadrilat- 
Paes tO avoid violating the contours. Nete iso that the 
Searse-tc-fine transition for the torsion function response 
is fairly smooth whereas the ‘transition for the strain 
energy density and shear stress refinement tends to produce 
distorted and elongated elements. This is aggravated py the 
fact that, unlike the torsion function, the shear stress and 
strain energy density are not monotonic over «he domain. 

The solution results obtained for each grid are 
presented in the upper half of Table IV. Reatlosteqlanc=, 
the results of the refinements are disappointing in ccmpar- 
ison +0 the parenthetic values obtained uSing a uniforn 
eed. While all three criteria produce improved accuracy for 
the maximum shear stress, the errors for the maximum torsion 
function value grow extremely larg: Recallizg the observa- 
tions made for the cne-dimensional study, it would follow 
that such behavior is prebably due to the unusually lazge 
elements near the center of the domain. The entries in the 
Mewer half of Table IV eflect the drasti improvement 
obtained by simply ae a faw additional degrees-of- 
freedom along those element edges which grew exceptionally 
long during the optimization process, thus increasing their 
polynomial order from one to two. Not only did this modifi- 
cation reduce the large errors for the maximum torsion 
function estimation, but modest improvements were also 


mp 
Memasned in the sstimatzons of the other resultants as well. 


Sis) 





"Wa psasIJ-Jo-saaibap 
JO Asdgqunu awes Jo ptih AeauTT wAOTTtun JOJ SSNT RA OTPSYRUSTCGy 


(90°Z-) 79°Z- (Z8°9-) 8Z°?C- (40°00) Lttr°Oo- 2L9 68 99— ss (pow) 


| 
| 
(99°L-) 96°L- (0ogs°S-) S9°0- (90°00) §62°0 LL 80 7ge.s«é(pou) =aas 
| abs 
(Sovt-) Gre_- =—s- (oe S-) hzet- = (90°0) 6Z°0- 6L  OLL 2a (pow) (fr 


(Z8°L_-) Un eh- (oZ°9-) LR°E- (90°0) 9°0l OL OOL ras aas 
(77°Z-) 18°C- (OZ°L-) On°Z- (s0°O) h°Z 9S 18 99 Jt 
(9z°L-) B8°L- (oo°9-) LSG*L- (90°00) 0°S EL ZOt 78 (fr 
— e9/L  Oo/xeUp | xeu A JOd PON weT_ sinoyuoy 
“ON “ON “ON 
~IOIATG DATRIOETOSY shezx,usorzed 
S2[NSeYy uot yNTOS bhutanozuoy 


AI ATES 


SS LS A A Se TS wee meee lt ee ce > om ewe Se Cm mm ee rm, oe ee ees ou ee oe ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ae SP SS eee eS 





Pcewadain, the selection of the optimum grid wouid depend 
predominantly on the optimization goal of the analysis, but 
one would likely agree that the strain energy density varia- 
tion along with some modification to restrict excessive 
element growth provides the superior refinement criterion. 
An additional word of caution is in order for the 
contouring techniques for grid optimization. Because the 
problem must be completely redefined from scratch after the 
initial analysis, the praprocessing cost can become enor- 
mous, e@specially if several cycles are employed *o obtain 
more precise contours as some authors suggest. Unless there 
is available an interactive automatic mesh generator based 
on this technique, such as the ons described in Reference 
(Ref. 11], contouring shouid be abandoned in favor of some 
more easily implemented dpi) SpALNizZ ation sechniques 


employing Similar refrinement criteria. 


2. Selective Refinement 





The Simplest way to avoid the problems encountered 
in the contouring techniques is to perform the initial anal- 
ysis on a reasonably coarse grid and then to selectively 
refine these elements over which the strain energy density 
varies the most. The critical concern then arises as to how 
Coarse the initial grid should be. If the preprocessor 
employs the necessary constraint conditions to permit the 
"father-tc-four-sons" element subdivision scheme direcely, 
or if hierarchical refinement is employed, then the initial 
grid should b@ just coarse enough to adéquately @ 
problem and to limit the number of refinement cycles neces- 
Sary. The latter becomes even less 9f a concern if iterat 
solvers are emploved. If, on the othar hand, coarse-tc-f 
transition schemes are used to implament the h-versicn or 
Chiy lcw pceiynomial order elements are available in the 


Pevers On, then the initial grid must be sufficiently fine 
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so asncet to restrict severely the amount of refinement 
which can be performed in any given cycle. Unfortunately, 
the conditions under which this investigation was conducted 


were those of the lat*er. 
a. The h-Version 


Selective refinement Dy the h-version was 
performed on both linear and quadratic element grids. For 
the linear case, the initial analysis was performed ona 
uniform grid »sf 55 elements, 72 anodes, and 50 degrees-of- 
freedcm. The initial quadratic analysis employed an eight 
element uniform grid of 37 nodes and 20 degrees-of-freedon. 
The reason for such a qreat disparity in the number of 
elements for the initial analyses 1S tnat subdividing a 
quadratic isment introduces many more degrees-of-freedom 
than the subdivision crt a linear element. These numbers were 
chosen to provide approximately the same number of degrees- 
of-freedcm for the optimum grid cof th2 final cycle for each 
case. The initial analysis is performed and those elements 
Over which the strain ¢nergy density variation is signifi- 
cantly greater beceme candidates for refinement. The refine- 
ment is performed by subdividing each candidate element into 
moue New ones by censtructing 2 soarse-to-fine transition 
zone of "buffer" elements around the czefined region. 
Successive analyses and selective refinements are repeated 
until the maximum element strain energy density variation is 
approximately that of the remainder of the grid. The process 
is further improved when the nodal values of the strain 
energy density are used to indicate the general direction in 
which the refinement is to proceed. This permits muitiple 
refinements in the same cycle, thereby reducing the number 
Pmmeycic= LTequired te arrive at ths optimal grid. fer «his 
problem, the Linear ne required two refinement cycles 


while the quadratic id required three cyclas. The 
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selective refinement process is depicted in Figures 6. 


6.10 tol the linear and Gladratac element ang Boole 
respectively. 
The solution results for each selective refine- 


ment cycle are presented in Table V. The most impressive 
ieee Ae 


maximum shear stress estimate for successive cycles. Whil 


observation to be made is the significant improvement 


(D 


4 


there is also modest imprevement in the accuracy of th: 
torque estimate for successive cycles, when compated to the 
estimate cbteined frem the uniform grid of the same number 
of degrees-of-freedom and polynomial order, the refinement 
Peelmate cl torque is slightly poorer. This is because addi- 
tional degrees-of-freedom are heing introduced in only a 
small region of the domain but the torque, and energy, are 
global quantities. The author has no satisfactory explana- 
tion why the estimate for the maxinum torsion function 
improves for successive refinements of the linear grid but 
not for the quadratic case. However, as has already been 
mentioned, this particular solution parameter appears very 
sensitive te such preblem variables as nodal placements and 
eiem=ent Shapes; hence, its behavior is difficult +o predict, 
even when the refinement 1s applied to regions remots fron 
the point where the maximum torsion function value occurs, 
as was the cas@ in these exampl2s. FOr “COnpitat zonal 
reasons, itis desirable to crestrict the number of refine- 
ment cycles +o a necessary minimum. In this ¢xampie, the 
Guadratic grid required an additional cycle over the linaar 
@crid but this is because it is necessary ‘+0 perform the 
initial quadratic analysis using far fewer degrees-of- 
freedom. Therefore, the early cycles of the quadratic anai- 


ysis actually represent comparatively smaller problems. 
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b. The v-Version 


Before centinuing to ths next optimization tech- 
Rique, it is worthwhile to take a quick look at selective 
refinement employing the p-version. Because the finite 
element code uSed in this investigation only provided for 
element crders of one and two, ths advantages of the method 
cannot be fully realized, but the affects of a singie cycle 
can be examined. 

Beginning with thrse unitorn ‘epee aalol= 
differing numbers of iinear e¢lements, the initial anai 
were performed. In each case, the elements over which the 
strain energy density varied the most were transformed 
4-noded Lagrangian elements into 8-noded SeTERd2 Day 
elements by the addition of midsidée nodéss. The elem 
grids are shown in Figure 6.11 and the asterisks 4d 
those elements for which the polynomial order was incre 
In this e¢éxamole, the number of 2lsments to be refin 
each case was chosen so as to achisve approximately the same 
number of degrees-of-freedom after 2 single cycle. 

The solution results are snown in Table VI. 
Significant improvements in the astimate of the maxinun 
shear stress were achieved for each case. An improvement in 
the estimated torque was also realized for all three casés, 
but was more noticeable when the number c fined elements 
Was larger. This is because gquadrtatic e2iéements are far 
superior to linear elements in tné modeling of integral 
guantities, as observed in Figure 6.4. Somew 
is the increased error in the 23stimate of the nmaxinua 
*crsion function value observed in two of she three refine- 
ments even though the elements in the vicinity of its occur- 
rence were not affected. Again, this is likely attributable 
to the unusual behavior patterns of this quantity already 


discussed. 
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In order to perform additional cycles of whe 
p-version, it would be necessary to alter the refinement 
criterion slightly. Because the element sizes do not change 
for successive cycles, the need for refinement would neces- 
sarily be based on strain energy density variation between 
nodes rather than over the elements. 

Selective refinement employing ‘the p-version is 
most efficientiy implemented hierarchically, in which case 
it acquires some attractive computational advantages. I+ is 
unfortunate that time did not permit further investigation 


here, but the need for future research is évident. 


3. Subdonain Isoiation 





The refinement criterion and initial precedures in 
mploying subdomain isolation are identical to those used in 
ements for 


selective refinement. After the candidate ei 
refinement are identified, they ara completeiy is 
the remainder of the domain and solved as sma 
Mecodems. The advantages of the technique are twotcld. By 


isolating the elemerts *c be refined the solution is not 


repeated in each cycle for those elements for which the 
initial analysis is assumed adequate. Furthermore, by elin:i- 
nating most of the i ncte-ntete 2eedom over the entir 


domain, the subsequent refinement of the isolated redisc 


> 


an 
be much greater than wouic otherwis2? be practical. 

As before, the technique was applied to both linear 
epeeduedratic uniform elemant grids. Those ¢lements of the 
initial analysis over which the strain energy density varia- 
tion was exceptionally large were isolated to comprise the 
subdomain in each case. There were three such ele 
f@emeimiicial linear grid and two fcr the quadratic grid. Each 
subdomain wes uniformly refined tc achieve approximate 
Same number of degrees-of-freedom as the initial analyses. 
Mem PLOCEss is devicted in Fiaqure 5.12 for «he linear grid 


aioe eogure 6.13 for the quadratic case 
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Fagure 6.12 Subdomain Isolation = Linear Elements. 
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Meagure 6.13 Subdomain Isolation ~ Quadratic Elements. 
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Ha DeEtommang Subdomain isolation <he gcverning 
equations remain the same, while oniy the domain and che 
Deundary conditions are altered. when che subiomain 
Beumdary has nodes common to the initial gzid, <7 one ae 
boundary values for those nodes are simply thea soiution 
values obtained from the initial analysis. The boundary 
Values arising from the introduction of new boundary nodes 
during the refinement process must be ie DY Neer CO" 
Maron Of the soluticn%esults of the initial analysis. one 
Of the options for an interpolation scheme iS Simply to use 
+he shape functions cf the unrefined elements. This may not 


Llements, Sc a aigher 


he 


Bemaesiztabple inthe case of si ieisic es 


RH 
es | 
ct 
He 
}-2 
— 


order interpolation may be employed. example 2 third 


(D 
6 
{i 
Qa) 
:. 


order Lagrangian polynomial was convenient for the linsar 


( 


case since there are four nodes from the nitial an or 
along the right-hand boundary of the subdomain. Sin 

are Only two such nodes on the upper boundary, it is neces- 
Sary to “borrow some adjacent nodal values from the 
discarded portion of the domain in the generaticn of new 
boundary values using higher order interpolation. 

The solution results voresented in Table VII are 
guite remarkable. ieee. Singmer cycle, cne somgtion accuracy 
for the maximum shear stress has increased by a full crder 
of magnitude. No other optimization technique examined in 
this investigation produced such inprovemencts. Note thas 
Maem nigher order polynomial interpolation fot the boundary 
values did improve the solution results for the linear case. 
One cf the disadvantages of this technique is that the 
refinement can produce no improvemen= in the Staha= sen. ot 
local quantities outside the subdomain. AS in this example, 
Eee eStimation of the Haximim =fOrSiOy Eunction value 
Obtained from the initial analysis must b2 accented as the 
Optimum Since it occurs outside ¢he supdomain. Furthermore, 


Since its value is predominantly affected by refinements in 
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PaonseOLentgh dedaaenes, te is doubtfuz that isolating 
new subdomain using the elements adjacant to its sroint of 
occurrence would noticeably improve its accuracy. However, 
Since the torque is a globally computed quantity, refinement 
weer 6Cilmmprove the accuracy of its contrzbutisn from the 
subdomain resulting in improvement of its accutacy overail 
as observed in Table VII. I* is this strictly lecal nature 
of the subdomain isolation technigue which restricts its 
applicability. But aif ‘the optimization goals are well 
Seemed and GE is wnderstood under which conditions and for 
what parameters it is effective, it can be an extremely 


powerful grid optimization techniqus. 


4, Mesh 3zading 


JG2 
(t4 


While mesh grading is nearly always performed on an 
Ha-oriorit basis, it may aiso be employed adaptively to 
Smoveade ae Sinple grid opt®mizazion technigue. <aAfter an 
initial solution has been obtained, the analysis may be 


repeated using varicus combinations of grading ratios in 


order to achieve a more uniform distribution of the element 
Strain anergy density variations. Here the grading ratio 
re=ers to the constant ratio of adjacent element lengths 


along a beundary of the domain to which grading is applied. 
mene ase several drawbacks to the technique, che first 
Beama <=hat @ good combination of g=z=ading ratios nay be 
difficult to obtain ina rceasonabi2 number of cycles. The 
other _ is that 1£ smalisrc elements are desired in 
mote tha one regicn the domain must b2 refined and 
constructed by subregions. 


O 
ct 


Uimesetnacely, =h2 “deomaineot “his problem isn 
well suited for mesh grading since it possesses no favorable 
geometric symmetry. Hence, the resulting element slongation 
and distortion would become severe for larger grading 


BieLOs. For simplicity, the nodal placements will be biased 
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oniy along the two domain boundaries adjacent to the point 
of maximum stress and the same grading ratio will be 
for each. This will result in small elements near the botto 

of the keyway and large elements along the periphery of the 
Shaft. While this is not the most deésiratle grid topology, 
it will produce a more uniform distribution of the element 
SED variations. 

The technique was applied to poth linear and quad- 
ratic element grids starting with 4a uniform mesh and succes- 
Sively increasing the grading ratio until the elements along 
the shaft periphery exhibited SED variations as large as 
those for the elements along the kayway. In both cases this 
condition occurred beyond the point where excessive slement 
elongation would be expected to produce numerical inaccura- 
cles. Gradéd meshes for selected grading ratios r, are shown 
in Figure 6.14 for linear elements and Figure 6.15 for quad- 
ratic elements. The solution results are presented in 
Taole VIII. As can be observed, “h2 maximum shear strsass 
estimate improves for each successive increase in ths 


the 


OP 


grading ratio. However, the cost of such improvement 
accompanying degradation in <¢*thn2 aStimate of the maxiaua 


torsion function value. This is to be @xpected since the two 


’ 


Maxima occur 3t difterent Locations in the domain and there- 
fore decreasing the size of the elements in the vicinity of 
one willl necessarily increas? tha element size near the 
other. Note that this degradation is not asarly so savers 
Meee the quadratic celenents, and the accuracy actully 
h 


because the higher order interpolation can better accomodate 


(p 


ct 


improves tor a low valus of Peace “uatz oes Tais is 


larger elements. for both Lineaz andi quadratic element grids 


the optimal accuracy in torque estimation cccurs for the 


mojerately graded meshes. 
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TABLE VIII 
Mesh Grading Solution Results 


| 
| 
2 
| 
| 
| 
ew 


Percentage Relative Error 


Gradin 
Ratio 2 Wye Tmax /G9 I / Gé 
Linear Element Grids ( 78 elements, 98 nodes, 
72 degress-or-freedon ) 
Pen, 0.060 -6.06 -1.77 
1.1 De siey ) =2.63 -1.56 
2 0.3389 =1.03 -1.77 
lees 0.679 -0.4%84 -2.20 


AO ate, le cine FEE sills quem SUES ge gs a i Loe ee "= GD eee Seg 


mai A aie a yy a a cits te i i i I i a yt in Sls I ieee ee erty 


Quadratic FBFiement Grids ( 28 elements, 107 nodes, 

76 degrees-of-freedom ) 

0 -0.0093 -5.26 ~0.0116 

| 0.0063 -1.85 ~0.0064 
Nes.2 0.0162 -0.606 -0.0091 
hee. -0.0246 =~. 4 43 -0.0179 | 
re ee a 


ites tkKely chaz Some of the error is attributable 
Memetne Numerical inaccuracies dua to element distortion. 
When applying a grading technique the analyst shculid seek an 
Gquitable balance between the refinement cricterion and the 


eae CODCLOgY. 
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Pinay eeoince sc ne Grading ratio is usual 
to the nodal separation rather than the elen 
leagths, it is advisable to reposition the midside 
quadratic elements so that they lie near the center of the 
element edges. This will generally improve the accuracy of 
ali the solution parameters, especially if the grading is 


sonewhat extreme. 


E. GUIDELINES FOR TWO-DIMENSIONAL GRID OPTINIZATION 


Moe order tO precvzds= some guidelines for obtaining 
optimal finite element solutions for two-dimensional prob- 
ems it is helpful t¢ compare the solution results obtained 
fom this problem employing the optimization techniques 
available. Such a comrarison 1s prasented in Table IX. The 
Moot Portion is for those techniques for which the initial 
analysis was performed using linear elements and the lower 
POe2 On WSing quadratic elements. Note that all of the grids 
employ approximately the same number of degrees-of-freecon, 


which was the chosen measure of analysis cost in this 


investigation. 
In making this ccmparison lt i185 important to understand 
just how significant a change in error actually is. If the 


convergence rate of a solution parameter is very siow, even 
emeeoMall reduction in the error may requir many more 
degreées-cf-~freedom. For this reason, the numbers in paren- 
theses have peen included by 3ach error to provide a rough 
approximaticn of the number of degrees-of-freedom required 
to Obtain similar accuracy using a uniform grid of the same 
number of degrees-of-freedom and ¢leaments or the same poly- 
nomial order as the initial analysis. In some cases, the 
analyses were not actually performed but instead the numbers 
in parentheses were obtained by extrapolation of the error 
versus egilpaaaaiee Gubve (Fig. 6.4), and iqnoring 


Bemme-O5f and til-coOndaitioning ef¢ects. 
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The tit sct sopservartion to be made from Table IX i 


}? 
Q WW 
ie 
te 
th 
} 
‘ 


while all of the optimization techniques produced si 
cant improvements in the accuracy of the maximum derivative 
quantity Tmax, the same cannot be said for the maximun 
soiLution quantity Wax and the integral quantity TT. One 
might even conclude that grid optimization is not cost 
efrective in the ccmputation of these values since the 
uniform grid provides estimates which are nearly as accu- 
meee, and in Some Cases better, than the optimum grid. This 
Benciusion would be correct if the solution maximum and its 
integral were the only resultants of interest in performing 
the analysis. Since the purpose of this study was to find an 
Seramum grid which produced acceptable errors for all the 
resultants, the Unt orm grid is clearly inadequate. 
Moreover, since in the majority of engineering problems it 
is the derivative cf the solution variable which is of 
primary interest, it deserves special consideration in 
making this compariscn. 


Furthermore, one might conclude? that the reason the 


“ 


Sof in che faxinum solution variable is ves for the 
optimal grid is because the s“=rain enargy densi variation 
criterion always concentrates the degrees-of-freedom in the 
vicinity of ths maximum derivative value, which in tnis case 
does net coincide with that of the maximum solution variable 
value. Howevct, Such a “Conelusion is incorrect and might 


erroneously lead one to attempt refinement where the maxinun 


f 


BeeuerOn Variable occurs in an effort t0 improve its esti- 
Mate. Other investigations have revealed that in nearly all 
cases the maximum accuracies for all thre 
tants are obtained by refinement in the O 
strain energy density variation [Ref. 11]. A O 
felivece Of Increasing error for the optimal grid 
element distortion which was encountered in LOT ewe 


a 
techniques, selective p-version refinement and subdomain 


4 





weoobeaciOn. Skeh @d2storet5sn can be avoided but would require 
more sophisticated refinement techniques than were availanble 
mor this investigaticn. 

reasonable choice for the optimum grid in Table IX 
would be cneé ‘for which all three values in parentheses are 
as large as reasonably possible taking into consideration 
the number of cycles required to provide such accuracy. 
Based on such a criterion, the author is partial to subdo- 
Main isolation for the solution of two-dimensional cvrobiems 
using linear elements, and selective refinement for finite 
element solutions using quadratic elements. Clearly, before 
a concrete recommendation couid be 2242 for a wide range of 
applications, many more problems would have «o be studied, 
but these two tecnnigues were fairly simple to implement for 
@astandard finite elamext code and the accelerated conver- 
gence of the solution resultants of interest was impressive. 
Conceivably, even greater solution accuracies mignt be 
Sioa. ned by using two or more of these techniques in 
combination. 

Here again the crucial element in selecting *he proper 
Seramaication strategy is the precise definition of the 
purpose for which the finite elament analysis is to be 
pertormed. [NeCeecetaas Seon stabte TX tend to support the 


following recommendations and conclusions: 


‘1) Regardless of the optimization strategy chosen, 
higher order elements are indispensable if high 
Becieen@s 68SOS 96 asegsal Sscluticn cueantities are 


(2) If the nn 


aximum derivative of th2 solution variable 
is of greatest concern, rose 
1 


(v 
YW) 


peeeeiy Local cCrerine= 
NGmesubacomain 259la2ion technigues car 
tional accuracy for a minimum number of 
e 


Les. 
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(3) 


If the maxinum solution variabl2 value occurs ata 
point in the demain removed from tne vicinity or the 
maximum derivative values, then its best eStimate 
will likely be obtained using a reasonably fine 
uniform grid and selectively subdividing elements in 
the Eeqtens <6 large strain energy density 


variations. 
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VII. CONCLUSION AND RECOMMENDATIONS 

The purpose of this paper has been to present an over- 
view of some readily employed finite element grid cptinmiza- 
tion methods and to demonstrate their effectiveness in the 
application t> some Simple problems. This work is by no 
Means all inclusive and the subject is still in its infancy. 
While there are many competing approaches to the problen, 
there is much more research to be done before any one 
becomes widely accepted as a standard analysis tocl. Because 
of the limited time and resources available, some of the 
more sophisticated rafinement criteria and techniques which 
have been developed have not been examined in any detail. 
Instead, the approach has been +9 2xamine those techniques 
which can be easily incorporated in a basic finite element 
coc2. However, it is likely that sone recently developed and 
rather elaborate self-adaptive grid optimization codes will 
soon be available. 

Also, this paper has not examined the Laportant clesses 
@repecblems in dynani and monlin2ear analysis. There is 
considerable ongding research in the extension of these 
techniques tO such probiéns, but the increased complexity is 
evident. For examole, in vibration analysis there is an 
optimum grid for each unique eigenvalus, but it is for these 
types of preblems that grid optimization is most promising. 

At the beginning of this paper it was stated that the 
go¢l of grid s9ptimization was to obtain maximum solution 
Seemrac’y £05 4 given analysis cost. Throughou«= «his paper it 
has been shown that, prior to successfully embarking upon 
Such a strategy, the underlying purpose of the analysis must 
be explicitly defined. Hopefully, it has been demonstrated 


iitee, rsd Optimizaticn is bv no means an mréealistic goal 
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Sre@miuset@dr NObe attractive than the non-adaptive practices 
wicely used today. 
The fcllowing are recommendations for future research 


topics: 


(1) Investigation of more sophisticated refinement 


ig 
criteria based on element residuals and reliabie 
er 


Tor estimates. 


(2) Investigation aye sl optimization techniques 


employing adaptive application of the p-version. 


(3) Implementation of a finite s2ieament preprocessor for 


performing hierarchical grid refinement. 


(4) Impiementation of a self-adaptive finite element 


code. 


(5) Application of grid optimization techniques to prob- 


lems of dynamic and nonlinear analysis. 
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FORMULATION OF THE ELASTIC CABLE PROBLES 


A. PROBLEM STATEMENT 


Consider a perfectly elastic cable initially stretched 
between two fixed rccints a distance 2L apart and under 


tension T. If the cable bears a distributed load per unit 





VOX) 


Figure A.1 Tensicn Cable Under Distributed Loading. 


length £(x) as shown in Figure che governing differen- 


Reale 
tial equation for the downward defaction v(x) is: 
Z 
V 
0 + £(x) = 0 Cin ae) 
SuEyect to the essential boundary condition: 


Y (xX = +L) = 0 (EC tease) 
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Té the distributed load is a supportive load provided bv 


paeenKles FCUNGation Cr Medulius k such -that 
f(x) = - k v(x) 


and if a concentrated load P is applied at the midspan, 


Equation A.1 becomes: 
q? 
aaa - k v =0 (Ean. As.3) 
subject to the natural boundary condition: 


#) ot > - £f2 (Eqn. A.4) 


Be PROBLEM SOLUTION 


The analytical sclution of th2 two-point boundary value 


preblem is: 


/2 


v (Xx) “3 Mcennewle cosh Ae) = Sinh \-x j (0-49, 4 1) 
( 


where A =a k/T : 


The finite element solution is obtained by the Galerkin 
formulaticn using the consistent rather than the lumped 


meploximation for she distributed loading. 


oy 


ae 





FORMULATION OF THE TAPERED BAR PROBLEM 


Aw. PROBLEM STATEMENT 


Consider 2 tapered bar of length L and constant nodulus 
of elasticity E fixed at one end. The cross-sectional area 
A(x) varies linearly from AL at th fixed end to Ay at the 
aE Let the bar be loaded axially by a concentrated tip 
Meea P, and 2 distributed ioad for which the intensity q (x) 


varies linearly from a at the fixed end to qi at the tip as 





Figure B.1 Tapered Bar with Applied Loads. 


omen n in Figure B.1. Die Geveraing  ditterential equation 
mom the axial displacement u(x) is: 


& [38 | =. G (0 <x $k) (Eqn. 8.1) 


oR 





subject to the essential boundary condition: 


a” (x6 = eC (eG Men se.) 


amd the natural boundary condition: 





du _ Pp 
dx x=” EAR (Eqns Sie) 
be PROBLEM SOLUTION 
Let a =1- AL 7A. and pe 8c 
& 
For F(x) = PB | q@{x) dx the solution is: 
xX 
x 
a8 F(x) 
ao) = a. A(x) a 
qt eL il P | ax 8 Qxe 
=O pis ty, 7 eS Beles es Bx 
ais { eee?) > (la) oC et loka) x re 
(OS x 


The finite element solution is obtained by the Galerkin 


qt 


formulation uSing «he consistent rather chan the Jumped 


Meee oximation for the distributed loading. 
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APPENDIX C 
FORMULATION OF THE TORSION PROBLEM 


A. PROBLEM STATEMENT 


Consider a solid circular Shaft of radius "a" and Shear 
modulus G, with a circular groove, or ksyway, of tadius b 
along a generator of the shaft with the cross-section shown 


Figure C.1 Cross-section of Shaft with Keyway. 


Memes ogure C.1. An applied torque I will produce an angle of 


Dd 
maaaet per wnit length 86. Ticme = iieworlin cONndlt.On” 2s 
aera 2f 2 torsicnal sttess function y» exists such that 


the shear stress comercnerts are: 


100 





The governing differential equation for the torsional stress 
Lunc tan U (Xp y) “1s; 


ae 2, 
afb , of (Eqn. C.1) 
ox? ay" 





Emogece tc the Dirichlet condition, = 0 on the boundary. 


Be. PROBLEM SOLUTION 


Mreewscolution Of Equation C.1 (Ref. 22: pp. 41-143] 1s: 


2 
=— ) - &(x2+y?) a (Same, Cea) 


W (x,y) = a(x-b? 





x2 +y2 


The maximum shear stress ocurs at point A (Fidq. 


js- “Ce lye sand 2S: 
Tmax = GO (2a - b) 


The applied torque computed from the area integral: 


r= 2c0 | ada 
A 


—= 


0 ; 
= | at-sa?e?-2b'y. + (2a%b+7ab?) sin | (Eqn C. 4) 


where @ = arcos (b/2a). 


Mae Strain energy per unit length of shafz= is: 


yt = ef dA = %TO (Sqn ceo) 
A 


oe 


— a 





Mace Velriataorak  soOrmulation of the finite element 
MweiscxeMacich 2S presented in detail in Chapter 6 of 


Reference 25. 
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